C========================================================================
C      
C                       DOUBLE PRECISION
C
C========================================================================
C  The following routine gets the regular solution of the following
C  differential equation. ECMN is the energy in Rydbergs.
C        2
C       d S   [ LN(LN+1)       V(X)          ]
C       --- - [ --------  + --------- - ECMN ] S =0
C         2   [   X*X           X            ]
C       dX      
C
C  Note that UCENTR(X) = V(X)/X and CNTFUG(X) = LN(LN+1)/X/X
      
      subroutine reglr8(ln,ecmn,eta,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,regcut,expcut,reg,jstart,jstop,phase)
c
      implicit real*8 (a-h, o-z)
c
      dimension jdouble(njdouble),ucentr(nx),reg(nx),gridx(nx),
     >   cntfug(nx)
      complex*16 phase

      jstart = jfirst18(ln,ecmn,gridx,nx,jdouble,njdouble,regcut)
      jstop  = jlast18(ecmn,gridx,nx,jdouble,njdouble,expcut)

      wnn=dsqrt(abs(ecmn))      
      if (jstart.eq.1) then
         s1 = 0d0
      else
         s1=appf18(ln,wnn,gridx(jstart-1))
      end if 
      s2=appf18(ln,wnn,gridx(jstart))

      call nmrvf8(ln,ecmn,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,s1,s2,reg,jstart,jstop)
      
      call match18(ln,eta,ecmn,jstart,jstop,gridx,nx,reg,phase)

      if (abs(reg(jstart)/regcut).lt.1e-4) then
         print*,'JSTART could be bigger. JSTART, REG(JSTART), LN, ETA,'
     >      //' ECMN, JSTOP',jstart,reg(jstart),ln,eta,ecmn,jstop
      end if 
      if (abs(reg(jstart)/regcut).gt.1e+2.and.jstart.gt.1) then
         print*,'JSTART should be smaller. JSTART, REG(JSTART), LN,'
     >      //' ETA, ECMN, JSTOP',jstart,reg(jstart),ln,eta,ecmn,jstop
      end if 
      return
      end
c""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
C  Riccati-Bessel function for small rho. Same as spherical Bessel * rho
C_______________________________________________________________________
      function appf18(ln,w,x)
      implicit real*8 (a-h, o-z)
      if (w.ne.0d0) then
         rho=w*x
         bessel=1d0
         do i=1,ln
            bessel=bessel*rho/dble(2*i+1)
         end do
         appf18=rho*bessel
      else
         appf18=x**(ln+1)/dble(2*ln+1)
      end if 
      return
      end
c=========================================================================
C  This routine matches the input REG solution to the apropriate
C  boundary conditions. It matches at two points of REG to the coulomb
C  functions F and H = G+iF. The two points are chosen where REG is not
C  too small.
C=========================================================================
      subroutine match18(l,eta,ecmn,jstart,jstop,gridx,nx,reg,phase)
      implicit real*8 (a-h, o-z)
      dimension reg(nx),gridx(nx)
      complex*16 phase,f1,f2,h1,h2
      complex*16 xx,eta1,zlmin,cfc(1),cgc(1),cfcp(1),cgcp(1),sig(1)
      data pi/3.1415927d0/
      
      mode1=12
      kfn=0
      if (ecmn.gt.0d0) then
         eta1=cmplx(dble(eta),0d0)
      else
         eta1=cmplx(0d0,-dble(eta))
      end if 
      zlmin=cmplx(dble(l),0d0)
      nl=1
      wnn=dsqrt(abs(ecmn))
      if (wnn.lt.1e-3) wnn=1d0
      dx=5d0*pi/4d0/wnn

      j=jstop
C  Find the smallest value of REG closest to REG(JSTOP)
      do while (abs(reg(j)).gt.abs(reg(j-1)).and.jstop-j.lt.20)
         j=j-1
      end do

C  Find the previous largest value of REG
      do while (abs(reg(j)).lt.abs(reg(j-1)).and.jstop-j.lt.40)
         j=j-1
      end do
C  First matching point
      j1=j
      if (j1.lt.jstart) j1=jstart+1
      ifai=0
 10   xx = cdsqrt(dcmplx(ecmn)) * dcmplx(gridx(j1))
      if (abs(ecmn).lt.1d-5) then
C  Here to evaluate the Coulomb G.F. at zero energy
         xx = dcmplx(dsqrt(8d0*gridx(j1)))
         kfn=2
         zlmin=2*l+1
      end if 
      call coulcc(xx,eta1,zlmin,nl,cfc,cgc,cfcp,cgcp,sig,mode1,kfn,ifai)
      if (ifai.ne.0.and.j1.ne.jstart+1) then
         ifai=0
         j1=jstart+1
         goto 10
      else if (ifai.ne.0) then
c$$$         print*,'IFAIL, ECMN in routine MATCH:',ifai,ecmn
         jstart=jstop+1
         return
      end if 
      if (abs(ecmn).gt.1e-5) then
         f1=cfc(1)
         h1=cgc(1)
      else if (eta.ne.0d0) then
         f1=cfc(1)*dsqrt(gridx(j1))*pi
         h1=cgc(1)*dsqrt(gridx(j1))*pi
      else
         f1=gridx(j1)**(l+1)/dble(2*l+1)
         h1=gridx(j1)**(-l)
      end if 

C  Choose the other matching point
      do while (gridx(j)+dx.gt.gridx(j1).and.j1-j.lt.20)
         j=j-1
      end do
      j2=j
      if (j2.lt.jstart) j2=j1-1
      xx = cdsqrt(dcmplx(ecmn)) * dcmplx(gridx(j2))
      if (abs(ecmn).lt.1e-5) then
C  Here to evaluate the Coulomb G.F. at zero energy
         xx=dcmplx(dsqrt(8d0*gridx(j2)))
         kfn=2
         zlmin=2*l+1
      end if 
      call coulcc(xx,eta1,zlmin,nl,cfc,cgc,cfcp,cgcp,sig,mode1,kfn,ifai)
      if (abs(ecmn).gt.1e-5) then
         f2=cfc(1)
         h2=cgc(1)
      else if (eta.ne.0d0) then
         f2=cfc(1)*dsqrt(gridx(j2))*pi
         h2=cgc(1)*dsqrt(gridx(j2))*pi
      else
         f2=gridx(j2)**(l+1)/dble(2*l+1)
         h2=gridx(j2)**(-l)
      end if 

C  Find real and imaginary parts of the phase
      phase=(f1*h2-f2*h1)/(reg(j1)*h2-reg(j2)*h1)
      const=abs(phase)
      phase=phase/const
      do j=jstart,jstop
         reg(j)=reg(j)*const
      end do
      return
      end

C========================================================================
C  This function returns the index of the first X for which the regular
C  solution is < REGSTART.
C========================================================================
      function jfirst18(l,ecmn,gridx,nx,jdouble,njdouble,regcut)
c
      implicit real*8 (a-h, o-z)
c
      dimension gridx(nx),jdouble(njdouble)
      w=sqrt(abs(ecmn))
      if (w.ne.0.0) then
         tmp=0.0
         do i=1,l
            tmp=tmp+log(float(2*i+1))
         end do
      else
         tmp=log(float(2*l+1))
         w=1.0
      end if
      xstart=exp((tmp+log(regcut))/(l+1))/w
      j=max(int(xstart*nx/gridx(nx)),1)
      do while (gridx(j).lt.xstart.and.j.lt.nx)
         j=j+1
      end do
C
C  make sure that we are not near a point where DX doubles
      do jd = 2, njdouble - 1
         if (abs(j-jdouble(jd)).lt.3) j = jdouble(jd) + 3
      end do       
      jfirst18=j
      return
      end
C=========================================================================
C  This function returns the index of the last X for
C  which EXP(-WNN*X) < EXPCUT
C========================================================================
      function jlast18(ecmn,gridx,nx,jdouble,njdouble,expcut)
c
      implicit real*8 (a-h, o-z)
c
      dimension gridx(nx),jdouble(njdouble)
      if (ecmn.ge.0.0) then
         jlast18 = nx
      else
         wnn = sqrt(-ecmn)
         xmax=-log(expcut)/wnn
         j=int(xmax*nx/gridx(nx))
         do while(gridx(j).gt.xmax.and.j.gt.1)
            j=j-1
         end do
         do while (gridx(j)**(1.0/wnn)*exp(-gridx(j)*wnn).gt.expcut
     >      .and.j.lt.nx)
            j=j+1
         end do 

C  make sure that we are not near a point where DX doubles
         do jd = 2, njdouble - 1
            if (abs(j-jdouble(jd)).lt.3) j = jdouble(jd) + 3
         end do       

         jlast18=j
      end if 
      return
      end
